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Problems consisting in finding the ground state of particles interacting with a given potential 
constrained to move on a particular geometry are surprisingly difficult. Explicit solutions have been 
found for small number of particles by the use of numerical methods in some particular cases such 
as particles on the sphere and to a much lesser extent on a torus. In this paper we propose a general 
solution to the problem in the opposite limit of very large number of particles M by expressing the 
energy as an expansion in M whose coefficients can be minimized by a geometrical ansatz. The 
solution is remarkably universal with respect to the geometry and the interaction potential. Explicit 
solutions for the sphere and the torus are provided. The paper concludes with several predictions 
that could be verified by further theoretical or numerical work. 
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I. INTRODUCTION 

Problems related to determining optimal particle dis- 
tributions under constraints are ubiquitous in the tradi- 
tional sciences and have been under intense scrutiny in 
the mathematical community 0, 0] • Two of the many in- 
teresting examples are the determination of ground states 
of M particles constrained to move on the sphere and in- 
teracting with a Coulomb potential, so called Thomson 
problem [3j, whose more direct representation are clas- 
sical electrons on helium bubbles Q and the crystalliza- 
tion of particles on spheres, relevant for understanding 
the structure of PMMA beads on oil/water droplets J5J- 
Similar problems on more general geometries such as the 
torus or negative curvature surfaces are also of great ex- 
perimental and theoretical interest 0, Q, @ ■ 

Theoretical investigations are surprisingly difficult. 
Extensive numerical results obtained in the Thomson 
problem, 0, 0, [n| , for example, show that the number 
of metastable states grow very fast with the total num- 
ber of particles, preventing a numerical solution to the 
problem even for a number of particles of a few hundred 
particles. For problems on the sphere, a few rigorous 
analytical results and conjectures on the energy of the 
ground state for large number of particles exist 0, H, 
but a description of the structure of these ground states, 
including practical tools on how to find them as well as its 
generalization to any given arbitrary geometry remains a 
completely open problem. 

Recently, it has been shown that elasticity theory 
H El El El (see also [H El) provides a power- 
ful framework to discuss best particle configurations on 
spheres, which can then be easily generalized to deal with 
any arbitrary geometry. Building on these results, we 
propose a general solution for the structure of ground 
states on arbitrary geometries in the limit of large num- 
ber of particles, and we construct the explicit solution for 
the case of a sphere and a torus. The solution is univer- 
sal, in the sense that it applies for short-ranged poten- 



tials and, in some situations, for long-ranged potentials 
as well. 




FIG. 1: Typical distortions due to geometric constraints on 
a perfect square lattice. If an additional row of atoms some- 
where in the middle of the crystal were added, variations of 
the lattice constant would become very small. 

In this paper, we will argue that the problem of finding 
the ground state of particles under constraints is equiv- 
alent to finding the particle distribution that is closer to 
a perfectly equilateral triangulation E3(t ne triangula- 
tion constructed from the actual distribution of particles 
via its Delaunay/Voronoi construction). The constraints 
that we consider in this paper are either geometrical (par- 
ticles are constrained on spheres, torus, etc..) or topo- 
logical (particles on a disk with the constraint that the 
total disclination charge is non-zero). If the particles are 
on a plane, the absolute ground state is a triangular lat- 
tice. The question is how the geometrical or topological 
constraints will modify this lattice. 

The effect of the constraints is to induce spatial vari- 
ations of the lattice constant, as it is sketched in fig- 
ureUfor a square lattice. This is the same situation that 
was first addressed by Frank^S] in the context of crystal 
growth, where he showed that in order to minimize the 
strains induced by the variable spacing it was necessary 
to add additional rows of atoms, that is, dislocations, 
which correct the spatial variations in the lattice con- 
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stant. Another critical result that we use in this paper is 
that the energy of particles interacting with a given po- 
tential can be regarded as an expression for large number 
of particles M 

whose leading coefficient contains 
a term, the C— function, which is always positive and en- 
codes the dependence on both geometry and topological 
defects. In this paper we will show how by adding dislo- 
cations it is possible to construct particle densities such 
that the C— coefficient becomes identically zero at lead- 
ing order in M (reaches its minimum value) thus pro- 
viding explicit ground states for very large number of 
particles. 

Although our arguments will be restricted to crystal- 
lization driven by energy, we believe that the geometric 
arguments used to construct ground states are general 
and apply to entropic driven crystallization such as for 
hard sphere potentials, since the entropic crystallization 
amounts to maximizing the area of the unit cell for a 
given packing fraction, and for larger packing fraction 
this leads to triangular lattices. 

The organization of the paper is as follows. In sec- 
tion [n] we review several results regarding the expansion 
of the energy for large number of particles. In section ITTT1 
we implement the solution outlined in fig. ^ for a gen- 
eral triangular wedge, and it is shown in section Hvl that 
the ansatz does minimize the energy. Explicit results for 
the sphere and the torus are provided in section We 
end with a summary of predictions that follow from this 
paper as well as with some conclusions in section IVTl 



II. ENERGY IN THE LIMIT OF LARGE 
NUMBER OF PARTICLES 

The energy of M particles interacting with a potential 
V(r) is given by 



E{M) = \Y,vm-?(j)) 



(i) 



where the sum runs over all M particles at positions 
f(i), (i — 1..M). For definiteness, we discuss the con- 
crete potential 



w \r\ s 
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but with minor modifications, the results can be made 
general to include any repulsive potentials. If the par- 
ticles are arranged in a configuration close to a trian- 
gular lattice, we write r(i) = R + u(i) + h(i), where 
R(i) — a(ne 1 +me 2 ) define the vertices of a triangular lat- 
tice of lattice constant a and primitive vectors e\, e2, and 
it, h are small quantities, in the sense that — << 1. The 
quantity u represents distortions tangent to the plane 
where h in the perpendicular direction. The energy Eq.^ 
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+ -jJ2ll (i,j)h(i)h(j)+O(u\u 2 h,h s )) (3) 

i,j 

where many higher order terms are neglected because of 
the assumed smallness of the displacements u, h. On a 
more rigorous basis, this step also requires the poten- 
tial to be short ranged (or s > 2), thus excluding the 
Coulomb potential, although for some geometries such 
as the sphere, we expect Eq. [21 to hold as well, as it will 
be shown. In the above expression, the contributions re- 
lated to the geometry and topological defects are entirely 
determined by the terms after the first. 

We now review how the energy can be regarded as 
an expansion in large number of particles. Detailed 
derivations have already been presented somewhere else 
0,0] , so we just recall the main results. To avoid exces- 
sive generality and keep the derivation simple, the results 
will be illustrated for the sphere. The first term in Eq.[3] 
gives the following explicit expression (for < s < 2) 
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where R is the sphere radius. The only modification for 
s > 2, is that the M 2 term, which arises from the long 
range nature of the potential, is absent. The actual val- 
ues for the function 6(s) maybe found in 0, Ej- The 
second contribution in Eq. is evaluated by retaining 
the leading term in an expansion in derivatives, leading 
to the familiar expression from elasticity theory 

2 2 
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where u a /3 is the strain tensor and A, \x are the Lame 
coefficients, whose explicit expression is ([HI) 



A* = 



r)(s) 



(6) 

for 0<s<2,^ = oo, but this is not important here (it 
becomes a constraint forcing the incompressibility of the 
crystal, but it has been shown that at leading order in 
M, the relevant elastic constant is the Young Modulus, 
which remains finite [l9|). It should be recalled that Eq.[3] 
involves terms which are not quadratic, so some higher 
order terms have been included. Combining Eq. [S] and 
Eq.^Jit follows 

d 2 r{ml + \{u a p?) = M 1+s/2 j^C(s, [u]) (7) 
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where the C— function is defined from 

G(s,[u])= (47r) l +s/2 J d 2 r{v{ S )ul + ^{u aP f). (8) 

The C— function has dimensions of area. Combining 
Eq. 0Eq. [5] and Eq.|U the energy Eq. becomes 
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Let us now discuss the dependence of the C— coefficient 
on the sphere radius R. For any configuration [it], we 
expect C(s, [u]) oc R 2 and this implies that the energy 
at order M 1+s / 2 is increased with respect to the planar 
value Eq. 0] If, however, one could find a configuration 
[u] such that its growth with R is linear at most, then it 
follows that 



C(s, u) oc Ra oc 



R 2 
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(10) 



where a is the lattice constant and in the last step, the 
_ 4-kR been used. In this case, for 
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M 



very large number of particles, the coefficient M 1+s / 2 in 
the energy expansion is given by the planar result and 
the configuration [u] becomes a minimum of the energy 
in the limit of large number of particles. 

We now analyze the approximation made in ignoring 
higher derivative terms in Eq. [5J Those terms will con- 
sist in higher derivatives of the strain tensor, which in 
turn will imply that the elastic constants, equivalent to 
the Lame coefficients, contribute to terms growing more 
slowly than M s / 2+1 , that is, the ignored terms do not af- 
fect the leading term in the expansion defined by Eq. 

In generalizing this expansion to other geometries there 
are several aspects to consider. First of all, the leading 
term for long-range potentials scales like M 2 . The gen- 
eral expression of the coefficient is 



E°= 6 4 [ d 2 xd 2 yp(x)-^- 
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where p{x) is the continuum density (The density of par- 
ticles at scales much larger than a). In a sphere, p is con- 
stant and given by p — , but in a general geometry, 
the density follows from minimizing the previous equa- 
tion under the constraint J d 2 xp(x) = M (For s = 1 this 
amounts to solving the Poisson equation for a fixed den- 
sity of charges). On a torus, for example, the resulting 
continuum density p(x) is not constant for s < 2. That 
is, the minimization of the leading coefficient for long- 
range potentials imposes density variations and invali- 
dates the form of the coefficient 0(M 1+S / 2 ) previously 
discussed. It is possible to generalize the expansion to 
include density variations but this will not be done here. 
On a general geometry, even for short-ranged potentials, 
the expansion Eq. 0] becomes slightly more complicated 



because there maybe more than one characteristic radius 
(In a torus, for example, there are two radii) and the coef- 
ficient C will contain a dependence on the dimensionless 
parameters that can be constructed from the geometry 
(for the torus this is the aspect ratio, the ratio of the two 
radii), but the property defining the ground state con- 
figuration is still given by Eq. 1101 The goal is now how 
to construct those configurations whose energy grows at 
most linearly with R. 



III. DISTRIBUTION OF DISLOCATIONS ON A 
TRIANGULAR WEDGE 

The problem consists now in obtaining the position 
and location of the dislocations needed to correct for the 
variable lattice constant, as discussed in the introduc- 
tion. This will be obtained from the following geometric 
argument. Let us consider a triangular patch of a sur- 
face like the one shown in figure [21 At one vertex (D 
in the figure), we place either a disclination of charge 
q = 1 (a vertex with five nearest-neighbors), q = (a 
regular six fold vertex) or q = — 1 (a vertex with seven 
nearest- neighbors) . 
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FIG. 2: (Color online) Triangle considered in the argument 



If we assume that a disk is formed out of identical tri- 
angular wedges, the angle 2a is given by 2a = ( for 



q = +1, 2a = =£■). We now discuss the necessary condi- 
tions so that the wedge may be triangulated with many 
equilateral triangles. For that matter, we now compute 
by how much the length of segment DV , which is given 
by r, differs from the length of the segment whose ori- 
gin is at V and forms 60 degrees with the segment DV 
(because we assume that the triangle whose vertex is at 
V is equilateral). The quantity to compute is l{r) — r, 
which from the purely geometric arguments outlined in 
fig. 121 for the geometry of a plane is given by 



l(r) — r = 2 sin( 



6(6-%) 



)h. 



(12) 



if l(r) — r = 0, which happens when qi = 0, the wedge can 
be tiled with all perfect equilateral triangles, but if the 
central disclination is non-zero, this is not possible. As 
explained in the introduction, we can fix this situation 



by including (or removing) additional rows of particles 
every time the equation 



l(r) 



±a 



(13) 



is satisfied, where a is the lattice constant. In physical 
terms we are adding a dislocation. In the planar case, the 
formula above implies that we add a grain boundary of 
equally spaced dislocations, where dislocations within the 



grain are separated a distance D 
suit is well known in metallurgy ! 



„ ■ , V r - This re- 

2sin{ 6(f ,'_ qi) ) 

I , where grain bound- 
aries of dislocations in the plane have been extensively 
investigated . 

In an arbitrary geometry, the function l(r) will also 
depend on the amount of Gaussian curvature enclosed 
within the triangular wedge, and this will lead to a dif- 
ferent type of grain boundaries. The function l(r) is now 
obtained from differential geometry, 



I(r) = / d^g ab (<p)V*(<p)V b (<p) 



*P\I ' 9ab{<p)V a {v)V b {v) (14) 



where g a b is the metric of the geometry, V a is the tangent 
vector of a geodesic, described by coordinates <p, which 
starts at point (r, — which is the equivalent of point 
V in figure [3 and forms an angle of 60 degrees with the 
direction defined by the geodesic defined by the radial 
distance r. Once the function l(r) is known, the location 
of the additional dislocations needed to ensure equilateral 
triangles will be obtained from the function l{r) — r from 
Eq. just as it was done for the planar case. In this 
paper, the distance l(r) is computed from a geodesic, but 
in general this is not necessarily the case, and examples 
will be given when discussing the torus. It should be 
pointed out that the predictions that follow from Eg. 1131 
and Eg. ll2l are equivalent to a similar formula provided in 
only in the plane, and differ on any other geometry. 
The simplest non-trivial example involving curvature is 
a spherical cap. We introduce a new parameter 9m that 
defines the angle subtended by the triangular wedge. We 
assume 



hi = L— 



R 



(15) 



where a is the lattice constant, L the total number of 
particles along the radial direction r and R the radius 
of the sphere. We consider a spherical cap both with a 
q = 1 disclination and without any disclination q = at 
its center, each case relevant for large and small aperture 
angle 9 respectively. 

The different functions l(r) — r are shown in figureOlfor 
the cases of a five- fold (q = 1) and (q = 0) disclination 
and compared to the equivalent results for the plane. The 
grain boundaries that follow from the function above and 
Eq. 1131 are shown in figure 01 for the particular situation 
L = 50. In the figure, squares represent seven-fold and 
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FIG. 3: (color online)i(r) — r on a sphere for both q = 1, 
and for a comparison, the results for the plane are also shown 



circles represent five-fold vertices, with dislocations being 
represented as a five-seven pair. For large aperture angles 
the last dislocation has rotated 180 degrees, and this is a 

reS 11 ^. nf thp fnnptinn l(r\ — r Vtppnmincr npcratwp at r I f? 

lai 



Spherical cap (L=50) 
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FIG. 4: (Color Online) Characteristic of grain boundaries of a 
spherical cap with different aperture angles with and without 
a disclination at the origin. Squares are seven-fold vertices 
and spheres are five-fold vertices. This plot corresponds to 
L/a = 50. 



Generalizations to any other geometry are now possi- 
ble. As an example, we discuss a wedge of a torus. The 
critical difference with the sphere is that the Gaussian 
curvature in a torus is not constant, and that implies 
that the function l(r) is different for the different wedges 
of a central defect. For a disk with a five- fold defect at its 
center, the five grain boundaries, which are identical for 
the sphere become now different. The second difference 
is that the line joining the dislocations is not straight, 
but curved, and depends on the torus aspect ratio 



Ri 
R 2 



(16) 
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where R\ and i?2 are the two torus radii. As an example 
the functions are shown for an aspect ratio r = 1.2. In 
this case. a.s shown in fiffiire^onlv three of the functions 



The energy Eq. El becomes 




Five-fold ( Curve I ) 
Six-fold (Curve 1) 
Five-fold (Curve 2) 
Six-fold (Curve 2) 
Five-fold (Curve 3) 



FIG. 5: (Color online) l(r) — r on a torus at r = 1.2 for both 
q — 1,0. If the central defect is oriented as shown in figure|H| 
only three functions are needed for q = 1 (fivefold). In this 
case R1/R2 = 1.2 and (3 = 0. 



IV. ENERGY OF OPTIMAL 
CONFIGURATIONS 

The previous argument provides the exact location 
where dislocations need to be added. Next, we compute 
the energy of these configurations. We will consider the 
energy for a large, yet finite number of part icles. The 
clastic energy Eq. can be discretized as [l4j 
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n, r \ + n- y (- - cos{9 df )) 2 . 
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where (b, c) runs over edges defined by nearest neigh- 
bors and angles. This energy has a very clear geometric 
interpretation as providing an energy cost for those tri- 
angles that either are formed of edges whose length is 
not the lattice constant, or whose angles are not exactly 
60°. The coefficients e, a are linearly related to the Lame 
coefficients |l ll l2l) . 

We now evaluate the energy function Eq.El m the case 
a = for the wedge discussed in the previous section. By 
construction, the length of the edges along the radial dis- 
tances r are exactly given by a, the lattice constant, and 
the contribution to the energy of these edges is zero, con- 
sistent with the total radial length of the wedge given by 
H = La, where a is the lattice constant and L the total 
number of edges. The main contribution to the energy 
then comes from the vertices along the direction defined 
by l(r). At radial distance r = ka there are n(k) of 
these vertices, and therefore, the average length of near- 
est neighbors along the radial distance ka is l(ka)/n(k). 



(l(ka) - n(k)a) 
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(18) 



The second term is proportional to the total number of 
dislocations Nri and takes into account that next to a 
dislocation the strains are significant. This energy in- 
cludes not only the core energy but also the (exponen- 
tially) small distortions arising from the grains. At this 
point it is illustrative to show the implications of the pre- 
vious formula for a simple and well known situation. We 
now compute the energy of a disk consisting of no defects 
other than a central disclination at the origin. In that 
case l(r) = br, where b is given from figure [21 and since 
no additional dislocations are added it is n(r = ka) = k. 
The above formula gives (with N4 = 0) 
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as expected, the energy of an isolated disclination grows 
quadratically with L 2 . Thus the C-coefficient defined by 
Eq. [8] will grow like H 2 , and lead to additional contribu- 
tions to Eq. [51 (here H = La plays the role of R in the 
sphere) . All these results are well known [14L l21j but it 
is of interest to show how they are recovered within the 
previous approach. It is now easy to show that the en- 
ergy of the configurations that satisfy Eq. ^| grow more 
slowly than R. We have, 
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where the last step follows because by construction 
l(r) — n(r)a never exceeds one. Within the same assump- 
tions, it is very simple to repeat the argument for the sec- 
ond term in Eq. 1171 and the same logarithmic behavior 
is found, so the statement holds for the entire range of 
elastic constants. As already mentioned, distortions near 
dislocations are not entirely negligible, and this leads to 
a term that grows linearly with L, . 



E x {e) = N d E c 



Ha = La 2 



(21) 



As an example of the previous considerations, we dis- 
cuss the spherical cap. We determine best distributions 
as a function of 6m, the subtended angle, and L the 
radial number of particles. As 9m — > 0, the spherical 
cap becomes a plane and the ground state should ap- 
proach the energy of a perfect planar lattice. As 6m is 
increased, Gaussian curvature effects become important, 
and a disclination at the top of the cap lowers the energy 
as shown in the inset of fig. El Upon including disloca- 
tions the energy decreases enormously and remains es- 
sentially independent with the subtended angle in both 
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FIG. 6: (Color Online) Energy of the spherical cap (L=50) 
with and without defects. The contribution proportional to 
the number of dislocations is not shown. 



The previous considerations do show that the geomet- 
rical ansatz does provide a minimum energy configura- 
tions, albeit a degenerate one. For the spherical cap, for 
example, the minimum energy can be achieved either by 
a plus disclination or without any disclinations at all, 
provided that the appropriate grain boundaries, as de- 
fined by Eq.^Hare used. The next question is therefore, 
which one of these minima is the actual ground state of 
the system. This involves considering the sub-leading co- 
efficients in the expansion. From Eg. I2UI this is given by 
the configuration with the lowest possible number of de- 
fects. This still does not completely solve the problem. 
One can consider several grain boundary on each trian- 
gular wedge where the separations of dislocations within 
the grain is larger, thus keeping the total number of de- 
fects constant. For the planar case, this question was 
investigated in [3] where it was concluded that grain 
boundaries with the smallest spacing within dislocations 
are favored. Further numerical investigations will hope- 
fully provide more evidence on this point. 



V. SOLUTIONS FOR THE SPHERE AND THE 
TORUS 

A. The Sphere 

There are several topological inequivalent triangula- 
tions of a sphere with icosahedral symmetry, and are la- 
belled by two integers (n,m). We just describe here in 
detail the solutions to the cases (n, n) and (n, 0), but so- 
lutions of the form (n, m) may be constructed along the 
same lines. The problem consists now in dividing the 
sphere into triangular wedges such that can be consis- 
tently joined back together after the additional disloca- 
tions necessary to relieve the geometric frustration have 
been included. 

The (n, n) configuration can be divided into 20 trian- 



gular wedges like the one shown in fig. (2 triangular 
wedges are shown). From the results described for spher- 
ical wedges, the dislocations follow the line AB, and the 
spacing is predicted from the function fig. [3] for an aper- 
ture angle 9 M = V = arcsin(2y / l/2 - l/(2v / 5)/v / 3) a 
37°. The fact that the angle ABC is 90 degrees and 
the angle ACB is 60 degrees, which can be checked by 
using formulas in spherical trigonometry, ensures that 
neighboring triangular wedges are perfectly joined and 
the complete consistency of the solution. 

The triangular wedges for the (n, 0) are shown in fig.0 
A natural patch is defined by DCO, but it should be 
noted that CO does not define a row of particles of 
the crystal. In this case, the patches extend beyond 
point C, and overlap slightly. The dislocations follow the 
lines DO and the spacing is predicted from the results 
from the function fig. [3] similarly as in the (n, n) case. 
Here, however, the consistency is a little more difficult 
to check, because in the region defined by such trian- 
gles such as COP, which spans aperture angles between 
(6 G ,0 V ), where 9 M = G = arccos(l/y(5))/2 a 32°, 
the different patches overlap. The consistency here is 
due to the sum of the Burgers vectors of the three grain 
boundaries that approach point O adding to zero, which 
ensures that all rows of particles added terminate within 
the patches. 




FIG. 7: (Color Online)Triangular wedges used to solve the 
(n, n) and the (n, 0) configurations. Dashed lines follow the 
directions of the grain boundaries. Filled circles represent + 
disclinations, also marked with D. 



B. The Torus 

The torus is depicted in fig. |H1 There are two radii 
curvature, i? 1 ,i? 2j and the aspect ratio is defined by 
r = jjf* > 1. As already discussed the only situation we 
discuss concerns short-ranged potentials. The Gaussian 
curvature of the torus depends only on the coordinate (p, 
and it is 



K = 



cos(ip) 



cos(yO) ' 



(22) 



Since Gaussian curvature attracts like sign disclinations 
@ , we assume here that the ground state of the torus con- 
tains 12 positive disclinations located along the geodesic 
of maximum curvature (outer curve) and 12 negative 
disclinations along the geodesic of minimum curvature 
(inner curve), and several grain boundaries of disloca- 
tions. Another possibility would consist of a ground state 
consisting of grain boundaries of dislocations only. It is 
very likely that for a thin torus r >> 1 disclinations 
may not be favored because nearest-neighbor disclina- 
tions are so far away, that radial grain boundaries may 
not efficiently screen the strains. 

If the ground state contains disclinations, the grain 
boundaries that appear following the five-fold disclina- 
tions have already been computed in figure [SJ There is 
a similar function for the seven-fold disclinations, which 
are located along the interior circle in the figure. The 
function predicting its spacing is given in figure 
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Seven-fold (Curve 3) 
Seven-fold (Curve 4) 
Seven-fold (Curve 5) 
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FIG. 9: (Color Online)Z(r) — r on a torus at r — 1. 
5 = 1. The different functions are shown. (R1/R2 
P = 0) 



2 for both 
= 1.2 and 



side of the torus. The particles are located in rings de- 
fined by the condition ip = constant, and therefore, only 
ip = 0, 7T are geodesies, which correspond to the two cir- 
cles in fig. Similarly as in the sphere, there are other 
torus triangulations that maybe constructed from twist- 
ing and rejoining the lattice either along the varphi or 
the 9 direction, but the discussion of the grain boundaries 
in this case is beyond the scope of this paper. 



VI. CONCLUSIONS 




FIG. 8: (Color Online)Representation of the torus, with 
the definition of the subtended angle. The two circles are 
geodesies of maximum and minimum curvature. A fivefold 
defect on the positive curvature of the torus has been shown. 



The situation were the ground state of the torus con- 
sists of grain boundaries of dislocations only is qualitative 
different, because in that case, the grain boundaries are 
not radial, but follow the directions defined by the az- 
imuthal angle, the curve AB in figure |S1 and a similar 
curve with opposite oriented dislocations on the inner 



We now summarize some of the main predictions that 
follow from the results presented in this paper that can 
be verified from either numerical simulations or experi- 
mental results: 

1. The ground state consists of grain boundaries of 
dislocations whose exact position and orientation 
is predicted from Eq. ^| Explicit examples have 
been provided for both the sphere and the torus. 

2. The Energy of the ground state tends to a universal 
value for very large number of particles, which is 
independent of the geometry and given by Eq. 0] 
(For long-range potentials, there is an additional 
term, which grows quadratically with the number 
of particles). 

3. For potentials V(r) — £7, the sub-leading correc- 
tions to the ground state are of order M( s+1 " 2 in 
the number of particles. 

4. The ground state is degenerate at leading order, 
as there are some free parameters characterizing 
the grain boundaries. This degeneracy is removed 
at sub-leading order by those configurations with 
the minimum number of dislocations and the lowest 
possible number of grain boundaries. 
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There are several ways to verify these predictions. The 
main difficulty with numerical minimizations is that it is 
very difficult to prepare an initial configuration that will 
relax to a previously selected distribution of defects. For 
the Thomson problem, howe ver, the ring-removal tech- 
nique developed by Toomre [22] seems to generate the 
type of grain boundaries that this paper predicts as the 
ground state. In fact, in [2j|, it was shown that these 
grain boundaries significantly lower the energy, bring- 
ing the icosadeltahedral configurations closer to the pla- 
nar limit Eq. It was further speculated that the true 
ground state could be achieved by successive applications 
of the technique, but this statement was not substanti- 
ated because the actual rings that needed to be removed 
were not known. The function in figure [21 together with 
Eq. E| does provide the location of the dislocations and 
therefore the actual rings that need to be removed, so 
the ring removal technique appears as a very promising 
practical tool to verify predictions 1 to 3 for the Thomson 
problem. It would also be of great interest to check how 
this predictions compare with ground states for short- 
ranged potentials, where the results of this paper are 
more rigorously justified. 

Evidence on the validity of statements 1-4 has also 
been provided from numerical and analytical results in 
|14|. where by using a discretized version of elasticity 
theory, it is shown that 1-4 are verified for a plane with 
the constraint of total disclination charge equal to ±1. 
The methods presented, however, are far general and it 
is expected that the extension of the results for any ge- 
ometry will soon follow. 

Statement 2) is in agreement with recent rigorous re- 
sults proven for potentials s > 2 ||. These results 
have greater generality, since they apply to dimensions 
other than 2. Preliminary results for the torus with short- 
ranged potentials have recently become available 0> Ull 
and show density variations for s < 2 and ground states 
for s > 2 with and without disclinations, but the present 
state of the simulations do not allow to draw a more 
quantitative analysis. Experimental evidence can also 
be used to prove the validity of statement 1. In [j| it 
was shown that PMMA beads assemble on an spherical 
oil- water interface forming spherical crystals, which can 
then be imaged by confocal microscopy. The current ex- 
periments show that next to five-fold defects additional 
dislocations arise, but there are only two of them, instead 
of the five predicted by this paper (Interestingly this is 
the number favored in the flat case 14] ). Furthermore, 
the dislocations within the grain show a constant spacing. 
This is possibly due to the fact that the total number of 
particles is still too small for the asymptotic results of 
this paper to apply. Future experiments with larger par- 
ticle aggregation numbers should settle this issue. 

An important point that has not been addressed in this 
paper is the critical value of M at which the asymptotic 



solution proposed will apply. A numerical verification of 
the predictions of this paper is currently under way for 
short-ranged potentials and preliminary results indicate 
that for the sphere, this number is of the order of five 
thousand particles at the very least. We hope to report 
more on this in the near future. 

In applying the formulas derived to a real situation, 
the dislocations must be located at points in the crystal, 
and accordingly, the spacing of the dislocations must al- 
ways be an integer. To apply the previous formulas, the 
location of the dislocations must therefore be rounded 
to the closer integer. The resulting orientations of the 
dislocations may also not be entirely consistent with the 
location of the crystallographic axis of the crystal, and 
this may require, for example, placing dislocations sepa- 
rated by an odd number of lattice constants 0] . Further 
numerical work will hopefully clarify this more technical 
points. 

There are a certain number of issues that this paper has 
not addressed. First of all, the ground states discussed 
for both the sphere and the torus rigorously apply only 
for some "magic" number of particles M. We expect that 
for large number of particles outside this magic numbers, 
the ground state solutions should not be that different 
from the ones proposed in this paper, since the range of 
M in between magic numbers is small compared with M 
itself. We hope that future numerical work will address 
this issue. 

In this paper, it has been assumed that the potential 
is isotropic. If the potential is not isotropic the ground 
state on a plane may not be a triangular lattice, and the 
present arguments need to be modified. For very large 
number of particles, it should be expected that locally, 
the triangulation will be very close to the flat case, so 
similarly as in the isotropic case, we expect that addi- 
tional dislocations will be required to fix the frustration 
induced by the geometric constraints. 

There has been recent interest in understanding sim- 
ilar problems as the one discussed in this paper where 
liquid crystalline order is discussed on a non-zero geom- 
etry. Examples include hexatic [2(| [53, UK , nematic (2t| , 
or smectic blue phases [3(3 . The results discussed in this 
paper, where the energy is regarded as an expansion in 
M and the defects correct for the geometric frustration 
are completely general and may be generalized to this 
cases as well. We hope to report more in the near future. 
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